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The existence and stability of localized patterns of criminal activity are studied for the reaction-diffusion model of urban 
crime that was introduced by Short et. al. [Math. Models. Meth. Appl. Sci., 18, Suppl. (2008), pp. 1249-1267]. Such 
patterns, characterized by the concentration of criminal activity in localized spatial regions, are referred to as hot-spot 
patterns and they occur in a parameter regime far from the Turing point associated with the bifurcation of spatially 
uniform solutions. Singular perturbation techniques are used to construct steady-state hot-spot patterns in one and two- 
dimensional spatial domains, and new types of nonlocal eigenvalue problems are derived that determine the stability 
of these hot-spot patterns to 0(1) time-scale instabilities. From an analysis of these nonlocal eigenvalue problems, a 
critical threshold K c is determined such that a pattern consisting of K hot-spots is unstable to a competition instability 
if K > K c . This instability, due to a positive real eigenvalue, triggers the collapse of some of the hot-spots in the 
pattern. Furthermore, in contrast to the well-known stability results for spike patterns of the Gierer-Meinhardt reaction- 
diffusion model, it is shown for the crime model that there is only a relatively narrow parameter range where oscillatory 
instabilities in the hot-spot amplitudes occur. Such an instability, due to a Hopf bifurcation, is studied explicitly for a 
single hot-spot in the shadow system limit, for which the diffusivity of criminals is asymptotically large. Finally, the 
parameter regime where localized hot-spots occur is compared with the parameter regime, studied in previous works, 
where Turing instabilities from a spatially uniform steady-state occur. 

Key words: singular perturbations, hot-spots, reaction-diffusion, crime, nonlocal eigenvalue problem, Hopf Bifurca- 
tion. 

1 Introduction 

Recently, Short et. al. [29, 30, 31] introduced an agent-based model of urban crime that takes into account repeat or 
near- repeat victimization. In dimensionless form, the continuum limit of this agent-based model is the two-component 
reaction-diffusion PDE system 

A t = e 2 AA - A + PA + a , d n A = , x e dfl , (1.1 o) 

rP t =DV- (vP-^-VA^j -PA + j-a, a;eO; d n P = , x € dfl , (1.16) 

where the positive constants e 2 , D, a, 7 and r, are all assumed to be spatially independent. In this model, P(x,t) 
represents the density of the criminals, A{x, t) represents the "attractiveness" of the environment to burglary or 
other criminal activity, and the chemotactic drift term — 2DV ■ {P-^jr) represents the tendency of criminals to move 
towards sites with a higher attractiveness. In addition, a is the baseline attractiveness, while (7 — a) /r represents 
the constant rate of re-introduction of criminals after a burglary. For further details on the model see [29]. 

In [29], the reaction-diffusion system (1.1) with chemotactic drift term was derived from a continuum limit of 
a lattice-based model. It was then analyzed using linear stability theory to determine a parameter range for the 
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Figure 1. Numerical solution of (1.1) at different times, for initial data close to a spatially homogeneous steady- 
state. Plots of A(x, t) are shown at the values of t indicated, (a) One dimensional domain with parameter values 
a = 1, 7 = 2, e = 0.02, r = 1, D = 1. Initial conditions are P(x, 0) = 1 — a/7 and A(a;,0) = 7(1 — 0.01 cos(67rx)). 
Turing instability leads to a formation of three hot-spots; one is annihilated almost immediately due to a fast-time 
instability, while the second hot-spot is annihilated after a long time, (b) D — 0.5 with all other parameters as in 
(a). Two hot-spots remain stable, (c) Numerical solution of (1.1) in a two-dimensional square of width 4. Parameters 
are a = 1, 7 = 2, e = 0.08, r = 1, D = 1. Initial conditions are P(x, 0) = 1 — a/7 and A(x, 0) = 7(1 + rand * 0.001) 
where rand generates a random number between and 1. 



existence of a Turing instability of the spatially uniform steady-state. A weakly nonlinear theory, based on a multi- 
scale expansion valid near the Turing bifurcation point, was developed in [30, 31] for (1.1) for both one and two- 
dimensional domains. This theoretical framework is very useful to explore the origins of various patterns that are 
observed in full numerical solutions of the model. However, the major drawback of a weakly nonlinear theory is 
that the parameters must be tuned near the bifurcation point of the Turing instability. When the parameters values 
are at an 0(1) distance from the bifurcation point, an instability of the spatially homogeneous steady-state often 
leads to patterns consisting of localized structures. Such localized patterns for the crime model (1.1), consisting of 
the concentration of criminal activity in localized spatial regions, are referred to as either hot-spot or spike-type 
patterns. A localized hot-spot solution, not amenable to an analytical description by a weakly nonlinear analysis, 
was observed in the full numerical solutions of [30] . 

As an illustration of localization behavior, in Fig. 1(a) we plot the numerical solution to (1.1) in the one-dimensional 
domain fi = [0, 1] with parameter values a = 1, 7 = 2, e = 0.02, r = 1, and D = 1. The initial conditions, consisting 
of a small mode-three perturbation of the spatially homogeneous steady-state A e — 7 and P e = (7 — a)/*) are first 
amplified due to linear instability. Shortly thereafter, nonlinear effects become significant and the solution quickly 
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becomes localized leading to the formation of three hot-spots, as shown at t w 14. Subsequently, one of the hot-spots 
appears to be unstable and is quickly annihilated. The remaining two hot-spots drift towards each other over a long 
time, until finally around t w 180, 000, another hot-spot is annihilated. The lone remaining hot-spot then drifts 
towards the center of the domain where it then remains. Next, in Fig. 1(b) we re-run the simulation when D is 
decreased to D = 0.5 with all other parameters the same as in Fig. 1(a). For this value of D, we observe that the 
final state consists of two hot-spots. Similar complex dynamics of hot-spots in a two-dimensional domain are shown 
in Fig. 1(c). 

It is the goal of this paper to give a detailed study of the existence and stability of steady-state localized hot-spot 
patterns for (1.1) in both one and two-dimensional domains in the singularly perturbed limit 

£ 2 «L>. (1.2) 

The assumption that e 2 -C D implies that the length-scale associated with the change in the attractiveness of 
potential burglary sites is much smaller than the length-scale over which criminals explore new territory to commit 
crime. In this limit, a singular perturbation methodology will be used to construct steady-state hot-spot solutions 
and to derive new nonlocal eigenvalue problems (NLEP's) governing the stability of these solutions. From an analysis 
of the spectrum of these NLEP's, explicit stability thresholds in terms of D and r for the initiation of 0(1) time-scale 
instabilities of these patterns are obtained. In a one-dimensional domain, an additional stability threshold on D for 
the initiation of slow translational instabilities of the hot-spot pattern is derived. Among other results, we will be 
able to explain both the fast and slow instabilities of the localized hot-spots patterns as observed in Fig. 1(a). 

In related contexts, there is now a rather large literature on the stability of spike-type patterns in two-component 
reaction-diffusion systems with no drift terms. The theory was first developed in a one-dimensional domain to analyze 
the stability of steady-state spike patterns for the Gierer-Meinhardt model (cf. [10, 3, 35, 37, 38, 34, 43, 46]) and, 
in a parallel development, the Gray-Scott model (cf. [4, 5, 15, 23, 24, 18, 1]). The stability theory for these two 
models was extended to two-dimensional domains in [39, 41, 40, 44, 43, 2]. Related studies for the Schnakcnburg 
model are given in [11, 36, 42]. The dynamics of quasi-equilibrium spike patterns is studied for one-dimensional 
domains in [9, 6, 7, 33, 22], and in a multi-dimensional context in [13, 14, 16, 2]. More recently, in [19] the stability 
of spikes was analyzed for a reaction-diffusion model of species segregation with cross-diffusion. A common feature 
in all of these studies, is that an analysis of the spectrum of various classes of NLEP's is central for determining 
the stability properties of localized patterns. A survey of NLEP theory is given in [46], and in, a broader context, a 
survey of phenomena and results for far-from-equilibrium patterns is given in [25]. 

In contrast, for reaction-diffusion systems with chemotactic drift terms, such as the crime model (1.1), there are 
only a few studies of the existence and stability of spike solutions. These previous studies have focused mainly on 
variants of the well-known Keller-Segel model (cf. [8, 12, 28, 32]). 

We now summarize and illustrate our main results. In §2.1 we construct a multi hot-spot steady state solution 
to (1.1) on a one-dimensional interval of length S. We refer to a symmetric hot-spot steady-state solution as one 
for which the hot-spots are equally spaced and, correspondingly, each hot-spot has the same amplitude. In §2.2 
asymmetric steady-state hot-spot solutions, characterized by unevenly spaced hot-spots, are shown to bifurcate from 
the symmetric branch of hot-spot solutions at a critical value of D. 

In §3 we study the stability of steady-state K- hot-spot solutions on an interval of length S when r = 0(1). A 
singular perturbation approach is used to derive a NLEP that determines the stability of these hot-spot patterns to 
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0(1) time-scale instabilities. In contrast to the NLEP's arising in the study of spike stability for the Gierer-Meinhardt 
model (cf. [35]), this NLEP is explicitly solvable. In this way, a critical threshold K c + is determined such that a 
pattern consisting of K hot-spots with K > 1 is unstable to a competition instability if and only if K > K c+ . This 
instability, which develops on an 0(1) time scale as e — > 0, is due to a positive real eigenvalue, and it triggers the 
collapse of some of the hot-spots in the pattern. This critical threshold K c+ > is the unique root of (see Principal 
Result 3.2 below) 



In addition, from the location of the bifurcation point associated with the birth of an asymmetric hot-spot equilibrium, 
a further threshold K c - is derived that predicts that a if-hot-spot steady-state with K > 1 is stable with respect to 
slow translational instabilities of the hot-spot locations if and only if K < K c _. This threshold is given explicitly by 
(see (3.20) below) 

K c _=( S -)p-^-«l!\ (1.4) 
V 2 J .Red v ' 



Since K c - < K c+ , the stability properties of a K- hot-spot steady-state solution with K > 1 and r = 0(1) are as 
follows: stability when K < K c -; stability with respect to 0(1) time-scale instabilities but unstable with respect to 
slow translation instabilities when K c - < K < K c+ ; a fast 0(1) time-scale instability dominates when K > K c+ . 

As an illustration of these results consider again Fig. 1(a). From the parameter values in the figure caption we 
compute from (1.3) and (1.4) that K c+ m 2.273 and K c - « 1.995. Therefore, we predict that the three hot-spots 
that form at t = 13.8 are unstable on an 0(1) time-scale. This is confirmed by the numerical results shown at times 
t = 25 and t = 30.8 in Fig. 1(a). We then predict from the threshold K c _ that the two-hot-spot solution will become 
unstable on a very long time interval. This is also confirmed by the full numerical solutions shown in Fig. 1(a). In 
contrast, if we decrease D to D — 0.5 as in Fig. 1(b) then we calculate from (1.3) and (1.4) that K c+ m 2.612 and 
K c _ ss 2.372. Our prediction is that the three hot-spot solution that emerges from initial data will be unstable on 
an 0(1) time-scale, but that a two-hot-spot steady-state will be stable. These predictions are again corroborated by 
the full numerical results. 

In §4 we examine oscillatory instabilities of the amplitudes of the hot-spots in terms of the bifurcation parameter r 
in (1.1). From an analysis of a new NLEP with two separate nonlocal terms, we show that an oscillatory instability of 
the hot-spot amplitudes as a result of a Hopf bifurcation is not possible on the regime r < 0(e _1 ). This non-existence 
result for a Hopf bifurcation is in contrast to the results obtained in [35] for the Gierer-Meinhardt model showing 
the existence of oscillatory instabilities of the spike amplitudes in a rather wide parameter regime. However, for the 
asymptotically larger range of r with r = 0(e~ 2 ), in §4.1 we study oscillatory instabilities of a single hot-spot in 
the simplified system corresponding to letting D — > oo in (1.1). In this shadow system limit, we show for a domain 
of length one that low frequency oscillations of the spot amplitude due to a Hopf bifurcation will occur when r > r c 
where 

r c ~ 0.039759(7- a) 3 a~ 2 £~ 2 • 

In §5 we extend our results to two dimensional domains. We first construct a quasi-equilibrium multi hot-spot 
pattern, and then derive an NLEP governing 0(1) time-scale instabilities of the spot pattern. As in the analyses of 
[39, 40, 41, 42, 43, 44] for the Gierer-Meinhardt and Gray-Scott models, our existence and stability theory for 
localized hot-spot solutions is accurate only to leading-order in powers of —1/ logs. In §5.1, we show from an analysis 
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of a certain NLEP problem that for r = 0(1) a if -hot-spot solution with K > 1 is unstable when K > K c where 

K c ~ ^" 1/3 l"l(7-a)^ 3 g - 4 /3 ( _ log£) V3 w (0 . 070 3 7) £> -i/3 a -2/3 (/y _ a) | | £ -4/3 ( _ log£) i/3 ; (L5) 

as e — »■ 0. Here to is the radially symmetric ground-state solution of Aw — w + w 3 = in R 2 and |f2| is the area 
of fi. As an example, consider the parameter values as in Fig. 1(c), for which |f2| = 16. Then, from (1.5) we get 
K c w 44.48. Starting with random initial conditions, we observe from Fig. 1(c) that at t = 9000 we have K = 7.5 
hot-spots, where we count boundary spots having weight 1/2 and corner spots having weight 1/4. Since K < K c , 
this is in agreement with the stability theory. 

Finally, in §5, we contrast results for Turing instabilities and Turing patterns with our results for localized hot-spots. 
We also propose a few open problems. 



2 Asymptotic Analysis of Steady-State Hot-Spot Solutions in 1-D 

In the 1-D interval x € [— M]i the reaction-diffusion system (1.1) is 

A t = s 2 A xx - A + PA + a (2.1a) 

rP t = D(p x -^-A^j -PA + ^-a, (2.16) 

with Neumann boundary conditions P x (±l,t) = A x (±l,t) = 0. Since P x — ^j-A x = (P/A 2 ) X A 2 , it is convenient to 
introduce the new variable V defined by 

V = P/A 2 , (2.2) 

so that (2.1) transforms to 

A t = e 2 A xx -A + VA 3 + a , (2.3 a) 

r(A 2 V) t =D(A 2 V x ) x -VA 3 +-/-a. (2.3 6) 

To motivate the e-dependent re-scaling of V that facilitates the analysis below, we suppose that D » I 2 and 
we integrate the steady-state of (2.3 6) over — I < x < I to obtain that V = c/ J_ A 3 dx, where c is some 0(1) 
constant as e — >• 0. Therefore, if A = (D(e~ p ) in the inner hot-spot region of spatial extent 0(e), we conclude that 
f_ l A 3 dx = 0(e 1 ~ 3p ), so that V = 0(e 3p_1 ). In addition, from the steady-state of (2.3 a), we conclude that in 
the inner region near a hot-spot centered at x = xo we must have —A + A 3 V = 0(A yy ), where y = (x — x )/e 
and A = 0(e~ p ). This implies that — 3p + (3p — 1) = —p, so that p = 1. Therefore, for D > ( 2 we conclude that 

V = 0(e 2 ) globally on — I < x < I, while A = 0(e -1 ^ in the inner region near a hot-spot. Finally, in the outer region 
we must have A = 0(1), so that from the steady-state of (2.3 6), we conclude that D (A 2 V X ) ~ a — 7 = 0(1). Since 

V = 0(e 2 ), this balance requires that D = 0(e~ 2 ). Since V — 0(e 2 ) globally, while A = 0(e _1 ) in the core of a 
hot-spot, we conclude that within a hot-spot of criminal activity the density P = VA 2 of criminals is 0(1). 

In summary, this simple scaling analysis motivates the introduction of new 0(1) variables v and D defined by 

V = e 2 v, D = D /e 2 . (2.4) 
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In terms of (2.4), (2.3) transforms to 

A t = s 2 A xx - A + e 2 vA 3 + a , -l<x<l; A x (±l,t)=0, (2.5 a) 

re 2 (A 2 v) t = D (A 2 v x ) x -e 2 vA 3 +7 -a, -I < x < I ; v x (±l,t) = 0. (2.5 6) 



2.1 A Single Steady-State Hot-Spot Solution 

We will now construct a steady-state hot-spot solution on the interval —l<x<l with a peak at the origin. In order to 
construct a K- hot-spot pattern on a domain of length S, with evenly spaced spots, we need only set I — S/ (2K) and 
perform a periodic extension of the results obtained below on the basic interval — I < x < I. As such, the fundamental 
problem considered below is to asymptotically construct a one- hot-spot steady-state solution on — I < x < I. 
In the inner region, near the center of the hot-spot at x = 0, we expand A and v as 

An 

A = — + A x H , v = v + evi + e 2 v 2 + e 3 v 3 H , y = x/s . (2.6) 

From (2.5 a) we obtain, in terms of y, that Aj(y) for j = 0, 1 satisfy 

Aq - An + v Al = , -00 < y < 00 , (2.7 a) 

A" - Ai + +3AlAivn = -a - viAl , -00 < y < 00 . (2-7 6) 

In contrast, from (2.5 b), we obtain that Vj for j = 0,1 satisfy 

(A 2 v' )' = 0, (A 2 v' 1 +2A A' 1 v' )' = 0, -00 < y < 00 . (2.8) 

In order to match to an outer solution, we require that Vn and V\ are bounded as \y\ — > 00. In this way, we then 
obtain that t>o and Vi must both be constants, independent of y. 

We look for a solution to (2.7) for which the hot-spot has a maximum at y — 0. The homoclinic solution to (2.7 a) 
with ^4q(0) = is written as 

Mv) = vo 1/2 My) , (2.9) 

where w is the unique solution to the ground-state problem 

w"-w + w 3 = 0, -oo<y<oo; io(0)>0, iu'(0) = 0; w^O as \y\ -> 00 , (2.10) 
given explicitly by w = \[2 sechy. Next, we decompose the solution A\ to (2.7 b) as 

Ai = a -Tr^w - 3awi , 

where w\{y) satisfies 

LqW\ = w'l — Wi + 3w 2 wi — w 2 , —00 < y < 00 , (2-11) 

with w[(0) — and w\ — > as \y\ — > 00. 

A key property of the operator L , which relies on the cubic exponent in (2.10), is the remarkable identity that 

L a w 2 = 3w 2 . (2.12) 

The proof of this identity is a straightforward manipulation of (2.10) and the operator L in (2.11). This property 
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plays an important role in an explicit analysis of the spectral problem in §3. Here this identity is used to provide an 
explicit solution to (2.11) in the form 

w\ = w 2 /3 . 

In this way, in the inner region the two-term expansion for A in terms of the unknown constants vo and vi is 

A(y)~e- 1 A Q (y)+A 1 (y) + -- - , A (y) = e-^'Mv) . Ai(y) = a (l - [w(y)} 2 ) - ^-w(y) . (2.13) 

In the outer region, defined for e <C |x| < I, we have that v = 0(1) and that A = 0(1). From (2.5), we obtain that 

A = a + o(l), v = h (x) +o(l), 

where from (2.5 b), ho(x) satishes 

h 0xx = C = ( " ~ f < , 0<|a;|<;; h 0x {±l) = , 

subject to the matching condition that h — > v as x — > The solution to this problem gives the outer expansion 

u~/i (ar) = ! (Z-|a;|) 2 -Z 2 + u , < \x\ < I. (2.14) 

Next, we must calculate the constants vo and v\ appearing in (2.13) and (2.14). We integrate (2.5 b) over —l<x<l 
and use v x — at x = ±1 to get 

e 2 j vA 3 dx = 2l(a - 7) . 

Since ^4 = 0(e _1 ) in the inner region, while A — 0(1) in the outer region, the dominant contribution to the integral 
in (2.1) arises from the inner region where x — 0(e). If we use the inner expansion A = s~ 1 Aq + A\ + o(l) from 
(2.13), and change variables to y = s~ x x, we obtain from (2.1) that 

/oo / />oo />oo \ 

A 3 dy + ehv A^Aidy + V! A 3 + 0(e 2 ) = 2l(a - 7) . (2.15) 
-00 V J —oo J —oo J 

In (2.15), we emphasize that the first two terms on the left-hand side arise solely from the inner expansion, whereas 
the 0(e 2 ) term would be obtained from both the inner and outer expansions. By equating coefficients of e in (2.15), 
we obtain that 

00 

2 



v Q = 2l(-f-a)/ A 3 dy, Vl = -3v A%A x dyj\ A 3 dy, 

J — oo J —oo J — oo 

Then, upon using (2.13) for A and A\, together with w = \/2sechy, w 3 dy — V2tt, and w 4 dy = 16/3, we 
readily derive from (2.1) that 

7T 2 ,3/2 ( JZo( w2 - w4 ) 4^3/2 2aTT 2 

vo = 7? , vi = 6v ' a -55 — = v ' a = -— — . (2-16) 

2P{a- 1 ) 2 ° ^ I Zo^ 3 dy J * ° Z 3 ( 7 ~a)3 

We summarize our result for a single steady-state hot-spot solution as follows: 
Principal Result 2.1: Let e — > ; and consider a one-hot-spot solution centered at the origin for (2.5) on the 
interval \x\ < I. Then, in the inner region y = x/e = 0(1), we have 

A(y) = — ^—+a [ l + ^^w-w 2 ] + o(l) , v ~ v + ev 1 H . (2.17) 



In addition, in the inner region, the leading-order steady-state criminal density P from (2.1) is P ~ w 2 . Here w = 
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Figure 2. Steady-state solution in one spatial dimension. Parameter values are Dq = l,e = 0.05, a = 1,7 = 2, and 
x e [0, 1]. (a) The solid line is the steady state solution A(x) of (2.5) computed by solving the boundary value problem 
numerically. The dashed line corresponds to the first-order composite approximation given by (2.19) (b) The solid 
line is the steady state solution for v(x). Note the "flat knee" region obtained from the full numerical solution in the 
inner region near the center of the hot-spot. The dashed line is the leading-order asymptotic result (2.18). 



A(0) (num) I A(0) (asyl) I A(0) (asy2) I «(0) (num) I v(0) (asyl) I v(0) (asy2) 



0.1 
0.05 
0.025 
0.0125 



6.281 
12.805 
25.628 
51.145 



6.366 
12.732 
25.465 
50.930 



6.003 
12.369 
25.101 
50.566 



3.5844 
4.1474 
4.4993 
4.7039 



4.935 
4.935 
4.935 
4.935 



2.961 
3.948 
4.441 
4.688 



Table 1. Comparison of numerical and asymptotic results for the amplitude A max = .4(0) and for v(Q) of a one- hot- 
spot solution on [—1, 1] with Do = 1, 7 = 1, and a = 2. The 1-term and 2-term asymptotic results for A max and v(0) 
are obtained from (2.17). 



w(y) = -\/2sech y is the homoclinic of (2.10), while vq and v\ are given in (2.16). In the outer region, O(e) < \x\ < I, 
then 



A ~ a + o(l); 



C 



((I - \x\ f - I 2 ) + V0 + 0(1) , 



(a - 7) 



< 0. 



(2.18) 



Note that to get a solution for A which is uniformly valid in both inner and outer region, we can combine the 
formulas (2.17) and (2.18). The resulting first-order composite solution is given explicitly by 



/2Z(7-a) 

\ 7T£ / 



a sech 



(!) 



+ a. 



(2.19) 



For a specific parameter set, a comparison of the full numerical steady-state solution of (2.5) with the composite 
asymptotic solution (2.19) is shown in Fig. 2. A comparison of numerical and asymptotic values for A(0) and v(0) at 
various e is shown in Table 1. From this table we note that the two-term asymptotic expansion for v(0) agrees very 
favorably with full numerical results. 
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2.2 Asymmetric Steady-State if-Hot-Spot Solutions 

In the limit e — > 0, we now construct an asymmetric steady-state if-hot-spot solution to (2.5) in the form of a sequence 
of hot-spots of different heights. This construction will be used to characterize the stability of symmetric steady-state 
-ftf-hot-spot solutions with respect to the small eigenvalues A = o(l) in the spectrum of the linearization. Since the 
asymmetric solution is shown to bifurcate from the symmetric branch, the point of the bifurcation corresponds to a 
zero eigenvalue crossing along the symmetric branch. To determine this bifurcation point, we compute v(l) for the 
one-hot-spot steady-state solution to (2.5) on |x| < I, where / > is a parameter. This canonical problem is shown 
to have two different solutions. A i^T-hot-spot asymmetric solution to (2.5) is then obtained by using translates of 
these two local solutions in such a way to ensure that the resulting solution is C 1 continuous. Since the details of 
the construction of the asymmetric solution is very similar to that in [36] for the Schnakenburg model, we will only 
give a brief outline of the analysis. 

The key quantity of interest is the critical value D^ K of D for which an asymmetric if-hot-spot solution branch 
bifurcates off of the symmetric branch. To this end, we first calculate from (2.14) that 

where the function B(z) on < z < oo is defined by 

B(z) = z 2 + 1/z 2 . (2.20 6) 

The function B(z) > in (2.20 6) has a unique global minimum point at z = z c = 1, and it satisfies B (z) < on 
[0,z c ) and B'{z) > on (z c ,oo). Therefore, given any z e (0,z c ), there exists a unique point z G (z c ,oo) such that 
B(z) = B(z). This shows that given any I, with l/q < z c = 1, there exists a unique [, with l/q = z > z c = 1, such 
that v(l) = v(l). 

We refer to solutions of length I and I as A-type and B-type hot-spots. Now consider the interval x € [a, b] with 
length S = b — a. To construct a JT-hot-spot steady-state solution to (2.5) on this interval with K\ > hot-spots 
of type A and Ki = K — K\ > hot-spots of type B, arranged in any order across the interval, we must solve the 
coupled system 2K\l + 2K 2 l = S and B(l/q) = B(l/q) for I ^ I. Such solution exists only if l/q < z c and l/q > z c 
with z c = 1. The bifurcation point corresponds to the minimum point where I = I = q. With D = D /s 2 , this yields 
that 

(Vir^f^f. (2 . 21 , 



\(j — a) 3 J \(j — a) 3 / 

At this value of the parameters, a steady-state if-hot-spot asymmetric solution branch bifurcates off of the symmetric 
if-hot-spot branch. This critical value of Dq determines the small eigenvalue stability threshold in the linearization 
of the symmetric i-C-hot-spot steady-state solution. For a symmetric configuration of K hot-spots on an interval of 
length S we have 2KI = S so that the critical value Do = D^ K , as defined by (2.21), can be written as 

A more detailed construction of the asymmetric solution branches parallels that done in [36] for the Schnakenburg 
model and is left to the reader. 
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3 The NLEP Stability of Steady-State 1-D Hot-Spot Patterns 

We now study the stability of the -ftT-hot-spot steady-state solution to (2.5) that was constructed in §2. The analysis 
for the "large" 0(1) eigenvalues in the spectrum of the linearization is done in several distinct steps. 

Firstly, we let A e , v e denote the one-hot-spot quasi-steady-state solution to (2.5) on the basic interval |x| < I, 
which was given in Principal Result 2.1. Upon introducing the perturbation 

A = A e + cj)e xt , v = v e + iPe xt , (3.1) 

we obtain from the linearization of (2.5) that 

e 2 4>xx -<P + 3e 2 v e A 2 e cp + e 3 Alip = \<j> , (3.2 a) 

D (eA 2 e ^ x + 2A e v ex <f>) x - 3e 2 A 2 e v e $ - e 3 ?/;A 3 = r\e 2 {eA 2 ^ + 2A e v e <p) . (3.2 6) 

We consider (3.2 a) and (3.2 b) on \x\ < I subject to the Floquet-type boundary conditions 

4>{i) = z<t>(-i), <f>'(!) = z<t>\-i), v(0 = ^H), v'(0 = ^'(-0. (3-2 c) 

where z is a complex parameter. 

For simplicity, in this section we will set r = in (3.2 b). The analysis of the possibility of Hopf bifurcations 
induced by taking r ^ is studied in §4. 

After formulating the NLEP associated with solving (3.2) for arbitrary z, we then must determine z so that we 
have the required NLEP problem for a -fT-hot-spot pattern on [—1, (2K — 1)1} with periodic boundary conditions. This 
is done by translating <p and ip from the interval [—1, 1] to the extended interval [—1, (2K — 1)1} in such a way that the 
extended <fi and ijj have continuous derivatives at x = 1,31, ... , (2K — 3)1. It follows that <p \(%K ~ 1)'] = z K (f>{— I), 
and hence to obtain periodic boundary conditions on an interval of length 2KI we require that z K — 1, so that 

Zj=e 2, ij /K^ j =0 ,...,A--l. (3.3) 

By using these values of zj in the NLEP problem associated with (3.2), we obtain the stability threshold of a K- 
hot-spot solution on a domain of length 2KI subject to periodic boundary conditions. The last step in the analysis 
is then to extract the stability thresholds for the corresponding Neumann problem from the thresholds for the 
periodic problem, and to choose / appropriately so that the Neumann problem is posed on [—1,1]. This is done 
below. This Floquet-based approach to determine the NLEP problem of a X-hot-spot steady-state solution for the 
Neumann problem has been used previously for reaction-diffusion systems exhibiting mesa patterns [22], for the 
Gierer-Meinhardt model [34], and for a cross-diffusion system [19]. 

We now implement the details of this calculation. The asymptotic analysis for e — > of (3.2) proceeds as follows. 
In the inner region with y = e~ x x, we use A e — 0(e^ 1 ) and v ex -C 1, to obtain from (3.2 b) that to leading order 
[to 2 ^!/]^ = in the inner region, where w is the homoclinic satisfying (2.10). To prevent exponential growth for ip 
as \y\ — > oo, we must take ip = ipo where ipo is a constant to be determined. Then, for (3.2 a) we look for a localized 
inner eigenfunction in the form 

0~$(2/)j y = e~ 1 x. 

Upon using the leading-order approximation A e ~ e~ 1 v Q ~ 1 ^ 2 w in (3.2 a), we obtain to leading order that $(y) satisfies 
Lo$ + ~T72 w3 ^° = A$ ' -oo<y<oo; L <P = $"-$ + 3w 2 $ , (3.4) 

V 
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with $ — > as \y\ — > oo. Here w is given in (2.16). 

In the outer region, away from the hot-spot centered at x = 0, we have A e ~ a and u = 0(1), so that (3.2 a) yields 



£ 3 a 3 



. A+1 *. (3.5) 

Then, from (3.2 &), together with A e <~ a, we obtain the outer approximation D [ea 2 ip x + 0(e 3 )] x = 0(e 3 ), which 
yields the leading-order outer problem 

ip xx = 0, < |^| < Z , (3.6) 
subject to the Floquet-type boundary conditions (3.2 c). The matching condition for the inner and outer represen- 
tations of tp is that linx^o ip{x) = i[)q, where ipo is the unknown constant required in the spectral problem (3.4). 
However, the problem for ip is not yet complete, as it must be supplemented by appropriate jump conditions for i[) x 
across x = 0. 

We now proceed to derive this jump condition. We first define an intermediate scale n satisfying e <C f] <C 1, and 
we integrate (3.2 b) over \x\ < n to get 

D {sA 2 e ip x + 2A e v ex <f>) \\ = f (Ze 2 Alv e 4> + e 3 ^A 3 ) dx . (3.7) 

We use the limiting behavior as x — > 0^ of the outer expansion to calculate the terms on the left hand-side of (3.7). 
From A e <~ a, (2.14) to calculate v ea: (0 ± ), and (3.5) to calculate ^(O^), we obtain that 

D (sA 2 e ^ x ) \\ ~ Doea 2 (V>*(0+) - ^(0")) = eD a 2 [V> x ] , (3.8 a) 

D (2A e v ex( f>) |% ~ 2D Q a [<P(Q+)v ex (0 + ) - 0(O"H X (O-)] = 4D a</»(0+) Uex (0+) ~ y^(0)( 7 - a)l . (3.8 6) 

Here we have defined [il>x\ = ipx(0 + ) — tp x (0~)- 

Next, since n ^> O(e), we can estimate the integrals on the right-hand side of (3.7) by their contributions from the 
inner approximation A e <~ e _1 w 1 ^ 2 w(y), ip <~ -00, <ft ~ ^(j/)) ancl w e ~ uo- in this way, we calculate 

I (3e 2 A 2 e v e <j) + e 3 ipAl) dx ~ 3e / w 2 ^dy+-^ j w 3 dy . (3.9) 

Upon substituting (3.8) and (3.9) into (3.7), we obtain the following jump condition for i[) x across x = 0: 

vf' 2 J w 3 dy-^-^-a)l)+3j w 2 <Pdy. (3.10) 

For the range A > — 1, we can neglect the negligible 0(e 2 ) term in the jump condition (3.10). In this way, the 
problem for the outer eigenfunction ip(x) is to solve 

V> xx = 0, 0<|x|<Z; V(0 = ^H). ip'(l) = z^'(-l), (3.11a) 

subject to the continuity condition 0>(O + ) = ^(0~) = V'o and the following jump condition across x = 0: 

/oo />oo 
w 3 dy, a 2 = 3 / w 2 $dy. (3.116) 
-00 J —00 

Upon calculating -0(0) = "00 from this problem, the NLEP is then obtained from (3.4). 

Upon solving (3.11) for ip(x), and evaluating the result at x = we get 



-3/2 , , / dv 



2,3/2 ! 



Ul^dy z 



(3.12) 
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Next, we use (2.16) for v and w 3 dy = V2ir to simplify VQ 3 ^ 2 ~ip . In addition, we use (3.3) to calculate 
(z - l) 2 

- '— = -2 + 2Rc(z) = -2 [1 - cos (2irj/K)} , j = 0, . . . , K - 1 . 



(3.13) 



Upon substituting these results into (3.4), we obtain the following NLEP for a if-hot-spot steady-state on a domain 
of length 2KI subject to periodic boundary conditions: 



Lo® - Xf 



J^w^dy 



D a 2 TT 2 



= A$, — oo < y < oo ; $^0, \y\ -> oo . 



1 + 



4/ 4 ( 7 - a) 3 



(1 - cos {2wj/K)) 



, j = 0,...,A--l. 



(3.14 a) 
(3.14 6) 



The final step in the analysis is extract the NLEP for the Neumann problem from the NLEP (3.14) for the periodic 
problem. More specifically, the stability thresholds for a if-hot-spot pattern with Neumann boundary conditions can 
be obtained from the corresponding thresholds for a 2if-hot-spot pattern with periodic boundary conditions on a 
domain of twice the length. To see this, suppose that is a Neumann cigcnfunction on the interval [0,a]. Extend 
it by an even reflection about the origin to the interval [—a, a]. Such an extension then satisfies periodic boundary 
conditions on [—a, a]. Alternatively, if (j>(x) is an eigenfunction with periodic boundary conditions at the edge of the 
interval [—a, a], then define 4>{x) — 4>{x) + <j){—x). Then, is a eigenfunction for the Neumann boundary problem on 
[0,a\. 

Therefore, to obtain the NLEP problem governing the stability of an steady-state X-hot-spot pattern on an interval 
of length S subject to Neumann boundary conditions, we simply replace cos(2-Kj/K) with cos(nj/K) in (3.14) and 
then set I = S/(2K) in the NLEP of (3.14). In this way, we obtain the following main result: 

Principal Result 3.1: Consider a K-hot-spot solution to (2.5) on an interval of length S subject to Neumann 
boundary conditions. For e — > 0, and t = 0(1), the stability of this solution with respect to the "large" eigenvalues 
A = 0(1) of the linearization is determined by the spectrum of the NLEP 

w 2 $dy 

LoQ-XjW*—^ — — =A$, -oo<y<oo; $^0, \y\ -> oo , (3.15 a) 



Xj 



1 + -T77 TvT o (1 - cos (irj/K)) 



4( 7 -a) 3 VS 1 



0,...,K-l, (3.15 6) 



where w(y) is the homoclinic solution satisfying w" — w + w 3 = 0. 

The stability threshold for Do is characterized by the largest possible value of Do for which the point spectrum of 
(3.15) satisfies Re(A) < for each j = 0, . . . , K — 1. In contrast to the typical NLEP problem associated with spike 
patterns in the Gierer-Meinhardt, Gray-Scott, and Schnakeneburg reaction-diffusion models studied in [4], [5], [10], 
[15], [35], and [36], [43], the point spectrum for the non-self-adjoint problem (3.15) is real, and can be determined 
analytically. This fact, as we now show, relies critically on the identity Low 2 = 3w 2 from (2.12). 
Lemma 3.2 : Consider the NLEP problem 



L ® - cw 3 w 2 §dy=:\$, -oo<j/<oo; $^0, \y\ -t oo , (3.16) 

J — oo 

for an arbitrary constant c corresponding to eigenf unctions for which f^w 2 ® dy ^ 0. Consider the range Re{\) > 
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-1. Then, on this range there is only one element in the point spectrum, and it is given explicitly by 



X 



/OO 
w 5 dy. (3.17) 
-oo 



To prove this we consider only the region Re(A) > — 1, where we can guarantee that |$| — > exponentially as 
\y\ — > oo. The continuous spectrum for (3.16) is A < — 1, with A real. To establish (3.17) we use Green's identity on 
w 2 and $, which is written as (w 2 L a § — $L w 2 ) dy — 0. Since L ® = cw 3 w 2 $ dy + A$ and L w 2 — 3w 2 , 
this identity reduces to 



[ w 2 ^dy(\-'S + c [ w 5 dy\=0, 

J — oo V J — oo J 



from which the result (3.17) follows. We remark that for the corresponding local eigenvalue problem Lo$ = v<&, it 
was proved in Proposition 5.6 of [3] that the point spectrum consists only of uq = 3 and the translation mode v\ = 
(with odd eigenfunction) , and that there are no other point spectra in — 1 < v < 0. When c = 0, we observe that 
(3.17) agrees with Vq. As a further remark, the result (3.17), when extrapolated into the region A < — 1, suggests 
that there is a critical value of c for which the discrete eigenvalue bifurcates out of the continuous spectrum into the 
region A > — 1 on the real axis. 

By applying Lemma 3.2 to the NLEP (3.15) we conclude that Re(A) < if and only if 

u> 3 dy 




,. x . =2, j = 0,...,K-l. (3.18) 
J-co w d V J 

In obtaining the last equality in (3.18) we calculated the integrals using w = -\/2secIiy. Since xo — 3 < Xi < Xi < 
■ • • Xk-\i a one- hot-spot solution is stable for all D , while the instability threshold for a multi hot-spot pattern is 
set by Xk-i- In this way, we obtain the following main stability result. 

Principal Result 3.2: Consider a K -hot-spot solution to (2.5) on an interval of length S with K > 1 subject to 
Neumann boundary conditions. For t — 0, and in the limit e — > 0, this solution is stable on an 0(1) time-scale 
provided that D < Dg K , where 

L = 2{ 1 - a f{S/2f 
n ° K - K^a 2 n 2 [l + cos(7T/K)}- [6AJ) 

In terms of the original diffusivity D, given by D = e~ 2 D , the stability threshold is Dj^ = £~ 2 Dg K when K > 1. 

Alternatively, a one-hot-spot solution is stable for all Dq > 0, provided that Dq is independent of e. 

Although we have not calculated the stability threshold for the small eigenvalues for which A — > as e — > in the 
spectrum of the linearization (3.2), we conjecture that this stability threshold is the same critical value D^ K of Do, 
given in (2.22), for which an asymmetric iiT-hot-spot steady-state branch bifurcates off of the symmetric iiT-hot-spot 
branch. This simple approach to calculate the small eigenvalue stability threshold, which avoids the lengthy matrix 
manipulations of [10], has been validated for the Gierer-Meinhardt, Gray-Scott, and Schnakcnburg reaction-diffusion 
models in [37], [36], and [18]. Since D$ K < D^ K for K > 2, we conclude that a symmetric Zf-hot-spot steady-state 
solution is stable with respect to both the large and the small eigenvalues only when D < Dq K . 

We make two remarks. Firstly, for the case of a single hot-spot where K — 1 we expect that the stability threshold 
for Do will be exponentially large in 1/e, and similar to that derived in [14] for the Gierer-Meinhardt model in the 
near-shadow limit. Secondly, we remark that the possibility of stabilizing multiple hot-spots for (2.1) is in direct 
contrast to the result obtained in the analysis of [32] of spike solutions for a Keller- Segel- type chemotaxis model with 



14 



T. Kolokolnikov, M. J. Ward, J. Wei 




Figure 3. Instabilities of a two-hot-spot steady-state solution induced by increasing D. Left: two hot-spots are stable 
with D — 15. Middle: two hot-spots exhibit a slow-time instability when D — 30. Right: there is a fast-time instability 
when D = 50. The parameter values are fixed at e = 0.07 a = 1, and 7 = 2, on the interval x E [— 1, 3]. The initial 
condition for the full numerical solution of (2.3) consists of two hot-spots that are perturbed slightly from the 
steady-state locations. 



a logarithmic sensitivity function for the drift term. For this chemotaxis problem of [32], only a one-spike solution 
can be stable. 



3.1 Numerical Results 

We now compare our stability predictions with results from full numerical solutions of (2.3). As derived above, under 
Neumann boundary conditions the thresholds on D for the stability of a symmetric if -hot-spot pattern on a domain 
of length 2KL are 

Dl ~ . D 0K = D$ K ( — 2 ) . (3.20) 

q 7T \1 + cos (n/K)J 

To numerically validate these thresholds, we choose e = 0.07, a= 1, 7 = 2, L = 1 and K — 2, so that we have an 
interval of length 5 = 4. For these parameters, our predicted stability thresholds are Dj| ~ 20.67 and s» 41.33, 
and our initial condition is a two-hot-spot solution with hot-spot locations slightly perturbed from their steady-state 
values. For our full numerical solutions of (2.3) we choose either D = 15, D — 30, or D = 50. Our stability theory 
predicts the following; the two hot-spots are stable when D = 15; the two hot-spots are unstable with respect to 
only the small eigenvalues when D = 30; the two hot-spots are unstable with respect to both the small and large 
eigenvalues when D = 50. The full numerical results shown in Fig. 3 confirm this prediction from the asymptotic 
theory. 



4 Hopf Bifurcation of if-Hot-Spot Steady-State Solutions 

In this section we study the spectrum of (3.2) for t > 0. This is done by first deriving an NLEP similar to (3.15). 
Since the analysis leading to the new NLEP is very similar to that in §3, we will only outline it here briefly. 

For t <C C(e 2 ), we get ip m ipo in the inner region x — 0(s), and hence (3.4) for the inner approximation $(?/) for 
<p remains valid. For r <C 0(e 2 ), we get to leading-order that ip xx = in the outer region < \x\ < I and so (3.6) 
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still holds. However, for t^O, the jump conditions (3.7)-(3.9) must be modified. In place of (3.7), we get 

rv rn 
D (eA 2 e 4> x + 2A e v ex <j>) |% = / {?>e 2 A 2 e v e cf> + e 3 ipAl) dx + e 2 r\ / (sA 2 e tp + 2A e v e cj>) dx . (4.1) 

J — r\ J —rj 

The left hand-side of (4.1) was estimated in (3.8), while the first two terms on the right-hand side of (4.1) were 
estimated in (3.9). We then use A e <~ e _1 i> 1 w(y), tp ~ tpo, 4> ~ &{y), an d v e <~ vq, to estimate the last term on 
the right hand-side of (4.1) as 



e 2 r\ / (eA 2 e ip + 2A e v e cf>) dx ~ s 2 t\ — w 2 dy + 2^ / w$ 

J —r) L ^0 J — oo J —oo 

Upon substituting (3.8), (3.9), and (4.2), into (4.1), we obtain that 



dy 



D ea 2 [V> x ] + 0(e 3 ) = e 



J — : 



:! r w Hd y +^- 

V 3/2 



f 

j — i 



w dy 



+ s 2 t\ 



/oo poo 
w 2 dy + 2^/vq / w$ dy 
-co J-oo 



(4.2) 



(4.3) 



which suggests the distinguished limit r = 0(e 1 ). Upon defining t = 0(1) by 

r = e-Vo , (4.4) 
(4.3) yields the jump condition (3.11 b) for ip across x — 0, where do, ai, and a 2 in (3.11 b) are to be replaced by 

/oo \ poo poo poo 

w 3 dy~— w 2 dy, a 2 = 3 / w 2 <Z> dy + 2y/Uo~T \ w<£> dy . (4.5) 
-oo ^0 J — oo J — oo J — oo 

With this modification of the coefficients in (3.11 b), the outer problem for tp is still (3.11). 
This problem is readily solved for tp(x), and we obtain that ip = -0(0) is given by 



' f x oo w 2 $dy s 

J^oo w3 dy 



2t XJvo 



X^oo w ^ dy 



D a 2 Tr 2 (z- l) 2 | 2t A 



8? 4 (7-a) 3 z Z(7-a) 



(4.6) 



J-oo w3 dy 

Finally, upon substituting (2.16) and (3.13) into (4.6), the NLEP problem for the Floquct problem on [—1, 1} follows 
from (3.4). As shown in §3, this problem allows us to readily determine the corresponding NLEP for the Neumann 
boundary condition problem on an interval of length S. The result is summarized as follows: 

Principal Result 4.1: Let r = 0(e^ 1 ) as e — > and consider a steady-state k-hot-spot solution on an interval of 
length S with Neumann boundary conditions. Define t c = 0(1) by t — e _1 S'(7 — a)T c / {AK). Then, the stability of a 
symmetric K-hot-spot steady-state solution is determined by the NLEP 

3 ( JZo w2 ® d y' 



J-oo w3 dy 

with $ — > as \y\ — ^ oo. Here, we have defined Xj, Xij> an d Pj by 



w$ dy = A$ , 



-oo < y < oo , 



(4.7 a) 



Xj 



Xij = (r c A)xj 



D a 2 7r 2 K 4 
4(7 - a) 3 



ft = 1 + 



-) (1 - cos (wj/K)) , j = 0,...,K-l. (4.76) 



This NLEP, with two separate nonlocal terms, is significantly different in form from the NLEP's derived for the 
Gierer-Meinhardt and Gray-Scott models studied in [3], [4], [5], [35], and [15]. 

Principal Result 4.2: There is no value of t c > for which the NLEP of (4-7) has a Hopf bifurcation. 

We note that there is a key step in the derivation of Principal Result 4.2 which relies on a numerical computation, 
see below. A completely computer-free derivation of this result is still an open problem. 
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Figure 4. A plot of the numerical result for p(/x), as obtained from (4.12). Note that is monotone increasing. 



Derivation of Principal Result 4.2: We use the notation Jhdy = h dy. Upon using Green's identity / w 2 LQdy 
J $Lw 2 dy and L w 2 — 3w 2 , together with (4.7 a), we obtain 



f w 5 dy 
J J w 3 dy 



J w 2 <S>dy-^- (^J w 5 dy^j (^J w<S> dy^j = X J w 2 $dy. 



Upon solving for J w 2 $dy in terms of / w&dy, and then substituting into (4.7 a), we get 



-f(X)w 3 J 



wQdy = A$ , /(A) = 



G X 3 ( Jw 5 dy 



(4.8 a) 



_Xij Xij(3- A) \Jw 3 dy 

We then simplify /(A) by using (4.7 6), together with J w 5 — (3/2) Jw 3 , to obtain 

1 _ 2^- 9 
/(A) _ r c A + ~ r c A(3 - A) ' 

Next, we observe that the eigenvalues A of the NLEP problem (4.8) are the roots of the transcendental equation 



(4.8 b) 



1 



/(A) 



— I w (L — A) 1 w 3 dy . 



(4.9) 



Upon recalling that L w = 2w 3 , we calculate 

J w (L - Xy 1 w 3 dy = i J w (L - A) -1 [(L - X)w + Xw] dy = i ^ + ^ ^ to (£ - A) -1 w dy . 

Substituting this result together with (4.8 b) and / w 2 = 4 into (4.9), we obtain that A is a root of 



- / w (L - A) wdy=—t- 

2 J T C X T c 



To determine whether a Hopf bifurcation is possible we set A = iXi in (4.10) and replace (L — A) -1 by (L — A) -1 = 
(Ll + Xf) (L + iXi). Then, upon comparing the real and imaginary parts in the resulting expression, we obtain 
that \i = X 2 and any Hopf bifurcation threshold r c must be the roots of the coupled system 

54 f ,~ N -i , 18 



w 



9 \-n 4/3,- 

io + Mj L wdy = + 



IV 



[ L l + m) w dy 



T c fl(9 + fl) 



(4.11) 
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Upon eliminating r c from (4.11), we obtain a transcendental equation solely for p = Xj > 0, 

20. fw ( L o + m) 1 L owdy 

p(p) = 3 - 20, - -g-n , where p(p) = 1 . (4.12) 

« J w{L 2 + p) wdy 

By using the identities Lou> = 2w 3 and L^w = (w + yw') /2, the limiting behavior for p(/x) is readily calculated as 

p(oo) = J „ " , = -; p(0) = - — 5-^- = -w 1.6461. 

fw*dy 3' J(L^w) 2 dy ^ + ^ 

Moreover, direct numerical computations of p{p) show that it is an increasing function of p, (see Fig. 4). On the 
other hand, for p > we have 3 — 20j — ^-p < 1 since Bj > 1. It follows that (4.12) cannot have any solution with 
p > 0. Consequently, there is no Hopf bifurcation on the parameter regime r = 0(e^ 1 ). ■ 

Such a non-existence result for Hopf bifurcations for the crime model when r = 0(e^ 1 ) is qualitatively very 
different than for the Gierer-Meinhardt and Gray-Scott models, analyzed in [35] and [15], where Hopf bifurcations 
occur in wide parameter regimes. 

4.1 A Hopf Bifurcation for the Shadow Limit 

Principal Result 4.2 has shown that there is no Hopf bifurcation for the regime t = C(e _1 ). However, a Hopf 
bifurcation can and does appear when r = 0(s~ 2 ). As will be shown below, in such a regime the amplitude of the 
hot-spot becomes oscillatory with an asymptotically large temporal period, due to an eigenvalue that is dominated, 
to leading-order in e, by its pure imaginary part. To illustrate this phenomenon, in this section we analytically derive 
the condition for a Hopf bifurcation of a single boundary spot on a domain of length one. To further simplify our 
computations, we will assume that Do in (2.5 b) is taken sufficiently large such that v(x, t) = v(t) can be approximated 
by a time-dependent constant. The limit D — > oo is called the shadow-limit (cf. [38]). Our main result is the following: 
Principal Result 4.3: Suppose that Do >• 1 and consider a half hot-spot of (2.5) located at the origin on the domain 
[0, 1], as constructed in Principal Result 2.1. Define to c by 

(24 - 7T 2 ) ( 7 - a) 3 (7 -a) 3 

t c= y QR 2 JU 2 ' = 0-039769 U ' . 4.13 
367H a 1 a z 

Let t = To/e 2 . Then, there is a Hopf bifurcation at r = t 0c . That is, the hot spot is stable for t < t 0c and is 

unstable for r > r 0c . Destabilization takes place via a Hopf bifurcation. More precisely, when r = 0(1), the related 

stability problem has an eigenvalue near the origin with the following asymptotic behavior as e — > 0: 



Numerical example. To illustrate Principal Result 4.3 we take 7 = 4, a = l,e = 0.05 and Do = 1000. Then, 
(4.13) yields r 0c w 1.07378. Now take r = 0.95 so that (4.14) yields the eigenvalue A w 0.4082i - 0.00529. We then 
expect the single hot-spot to be stable, although it will exhibit long transient oscillations. From the eigenvalue, we 
can estimate the period of the oscillation to be P — 2 ^ g2 m 15.39. This agrees with full numerical solutions of (4.15) 
as shown in Fig. 5(a). 

Next, we increase To to 1.15, while keeping the other parameters the same. In this case, tq > t 0c so that the 
hot-spot is unstable in the limit e — > 0. However, to = 1.15 is very close to the threshold value, and with e = 0.05, 
we expect even longer transients with the final state still unclear at t = 300. This behavior is shown in Fig. 5(b). 

Finally, as shown in Fig. 5(c), when we increase tq to 1.35, we clearly observe oscillations of an increasing amplitude. 
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Figure 5. maxu versus t with e = 0.05, a = 1,7 = 4, t = tq/e 2 with ro as given in the figure, (a) To < t c = 1.074; 
damping is observed, (b) To > r c but ro is very close to r c ; eventual fate of the oscillation is unclear, (c) ro > r c ; 
oscillations of increasing amplitude are observed. 



Derivation of Principal Result 4.3: We begin by re-writing (2.5) as a shadow system. For convenience, we also 
rescale A as A = u/e. In terms of this scaling, u — 0(1) in the interior of the hot-spot. Expanding v in powers of 
Dq 1 we then obtain to leading order that v(x,t) ~ v(x). We then integrate (2.5 b) and use the no-flux boundary 
conditions to obtain the following shadow-limit system oni£ [0, 1]: 



u t = e 2 u xx — u + vu 3 + sot . 



where we have defined \i by 



t I v / u dx } = u vi u 

Jo Jt £ Jo 



H = 7 — a . 



3 dx: 



u x (0,t) = u x (l,t) = 0, (4.15 a) 



(4.15 6) 



The shadow problem (4.15) is the starting point of our analysis. The corresponding steady-state system is 



= e 2 u xx — u + vu 3 + ea , fi = —v I W <Ir : 



u x (0,f) =u x {l,t) =0. 



(4.16) 



As shown below, in order to analyze the Hopf bifurcation it is necessary to construct the steady-state solution to 
two orders in e. To do so, we let y — x/e and we expand 



U = Uq + £U\ . 



v = v + evi + 



(4.17) 



— 1/2 

By substituting this expansion into (4.16), and equating powers of e, we obtain that uq = v n w, where w(y) is the 
positive homoclinic solution of w yy — w + w 3 = 0. In addition, u\ satisfies 

LqU\ = —a — w 3 viv Q 3 ^ 2 , (4-18) 

where the operator Lq is defined by Lqu = u yy — u + 3w 2 u. This operator has several key readily derivable identities, 

L (l) = -l + 3w 2 , L w = 2w 3 , L w 2 =3w 2 , 

which allows us to determine the solution u\ to (4.18) as 

(4.19) 



2 ^1 

u\ = a — aw ttttW ■ 



2v, 



3/2 
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Next, to determine vo and v\ we must calculate the integral in (4.16) for e <C 1. This yields that 

li= - I vu 3 dx <~ / (vq + evi) (uq + SuqUis) dy ~ v^ 1 ^ 2 / w 3 dy + e (^3uqUiv + viUq) dy + 0(e 2 ) . 
£ Jo Jo Jo Jo 

By equating coefficients of e, we get that Vq and t>i satisfy 

/>oo />oo />oo 

w _1/ ' 2 / w 3 dy = fj, , / 3uiw 2 + ^i?; ~ 3/ ' 2 / w 3 dy = . 
Jo Jo Jo 

Upon using the solution (4.19) for u\, the unknown «i can be determined in terms of a quadrature as 

1-5 J °°w 3 d y • 
The integrals defining vo and «i are then calculated explicitly by using 

/oo />oo />OC />OG -j^g 

w 2 dy = 4, / w 3 dy = w dy = V2ir , / w 4 dy = — , 
- oo J —oo J — oo J — oo 

which yields the explicit formulae 

7T 2 

u = yM" 2 , «i = -2a^ 2 M - 3 . (4.20) 
Next, we study the stability of this solution. For convenience, we extend the problem to the interval [—1, 1] by 
even reflection. We linearize (4.15) around the steady-state solution to obtain the eigenvalue problem 

,i/e _ 1 ,i/e 

= e 2 4>" -</>+ 2,u 2 v(j) + u 3 tp ; tA / (2<jmv + u 2 ip) dy = — / (2>u 2 v(j) + u 3 ip) dy , 

J-i/e £ J-i/e 

where the constant tp denotes the perturbation in v. Upon solving for ip we obtain 

rAe f^_{ £ /r 2(j)uv dy + [^{%3u 2 v4>dy 

e 2 c/)" -<t> + 3u 2 v(/) + u 3 tlj = X(f>; if, = {j§ . (4.21) 

T\eJ_ 1 ^ £ u 2 dy + J_ 1 ^ £ u 3 dy 

To motivate the analysis below, we first suppose that rAe 3> 1. Then, by using u ~ w/^/vq and v ~ vq , (4.21) 
reduces to leading order to the NLEP 

L <p-2w 3 l^ = \<f>. 
J w 

Here and below, J f denotes fdy. This problem has a zero eigenvalue corresponding to the eigenfunction <j> = w. 
All other discrete eigenvalues satisfy Re(A) < (cf. [46]). Therefore, the critical eigenvalue will be a perturbation of 
the zero eigenvalue. 

A posteriori computations shows that the correct anzatz is in fact 

A = e 1/2 A +eAi + ... ; t = t £~ 2 . 

The analysis below shows that Ao is purely imaginary, and hence determines the frequency of the oscillation, but 
not its stability. Therefore, a two-term expansion in A must be obtained in order to determine the stability of the 
oscillations. As such, we must expand all quantities in the shadow problem up to 0(e). The delicate part in the 
calculation is to note that 



where the last integral is in fact 0(e) as a result of 

,i/e nl/e 



f ' 2 f 2 f 2 / ' 2 

/ u = u + e 2u ui +e u t , 

J-i/e J J J-l/e 

as a result of 
f i/e rife 

/ u\=e 2 / (a + ...) 2 =2a 2 e + ... 
J -He J -lie 
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Thus, this term "jumps" an order and is comparable in magnitude to e J 2u u i- The remaining part of the analysis 
is more straightforward. We let t = tqe~ 2 and expand ip in (4.21) as 

rXe f 2(f)uv + f 3u 2 vcj) , i /o , 
^ = \ 1 2 Jr 3 = + £ 1/2 V-i/2 + eV>i , 

so that the eigenvalue problem for from (4.21) becomes 

(e 1/2 A + e\i)4> = Lq4> + (3ulvi + 6u uiv )4>e + uf, (ip + e 1/2 ip 1 / 2 + sipij + 3uluiip e 

= L^ + e 1/2 L 2 ^ + eL 3 <j). (4.22) 

After tedious but straightforward computations, the three operators in (4.22) are given by 

f web q 

L x 4> = L <j) - 2^-^w 3 ; 4.23 a 

J w 

L 2 <j) = VV2"o = (co j 'w<j> + a J w 2 cj^j w 3 ; (4.23 6) 

L 3 <f> = (3uq«i + 6u uiv )cf) + 3ului~ip + ufai 

= (c 2 w + c 3 w 2 + c 4 w 3 ) (j) + (c^w 2 + CqW 3 + c 7 w 4 ) J wef) + c 8 w 3 J (j) + c g w 3 J w 2 (j) , (4.23 c) 

in terms of the coefficients cq, . . . , eg defined by 



jt* 3^V2 3ttO!\/2 3V2 



co = -. — r- ; c\ = —- — ; c 2 = ; c 3 — ; c 4 = -c 2 ; c 5 = - 



7ra 



4toAo 47rr Ao /U 4^ 



"0 - 1 ^2 o 2x2 +~^~ ' C 7 = _C 5 ! c 8 = Z 5 Cg = — h -j- + j-j . (4.24) 

Finally, we expand in (4.22) as 

<j> = w + e 1/2 4>\ + s4>2 , 
and equate powers of e in (4.22). This yields the following problems for and <f) 2 : 

X w = Li4>i + L 2 w, (4.25) 

Xiw + Xo4>i — Li4> 2 + L 2 (j)i + L 3 w . (4.26) 

To determine Ao and <j>\ we must formulate the appropriate solvability condition based on the adjoint operator L\ 
of Li defined by 

L* l( p = L 0( p-2^ r ^- W . (4.27) 

Since L\ admits a zero eigenvalue of multiplicity one, then so does L\. In fact, w* defined by 

w* = (yw y + w)/2, (4.28) 

is the unique element in the kernel of L\, i.e. L\w* — 0, owing to the following two readily derived identities: 

f w 3 w* 



L w*= Wl 2 J , . = 1. (4.29) 



•j 3 w* 

Next, we impose a solvability condition on (4.25) in the usual way. We multiply (4.25) by w* and integrate by 
parts to derive that 
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where we used the integral identity in (4.29) together with J w*w = J w 2 /4. From the formulae for the coefficients 
Co and ci in (4.24), (4.30) determines Ao as 

\ Q = ±i.[^-. (4.31) 
V r o 

Since Ao is purely imaginary, the next order term Ai needs to be computed to determine stability. 
The problem (4.25) for <f>i can be written by using (4.23 a) and (4.23 b) as 

L 4>i = X w - ^c J w 2 + ci J w 3 + 2 y^r) w3 ■ 

Since L w* = w and L n w — 2w 3 , we can write <f>\ as 

0i = A w* - ^ where d = c f w 2 + c x f w 3 + 2 ^ ™ < t 1 . (4.32) 
2 J J J w 2 

Upon substituting (f>i into the definition of d, we can then solve for d by using (4.30) for Ao to get 
2d = c J w 2 + ci J w 3 + 2A f™™ 2 = c J w 2 + ci j w 3 + 2 (c a J w 2 + c x J w 3 ^ ^ 
Finally, by using J w*w = J w 2 /4, the expression above simplifies to 

d = Co J w 2 + ci J w 3 



,*„„3 



so that 0i is given explicitly from (4.32) as 



i = c / w + ci ; w 



(2w* - |) . (4.33) 



With (f>\ explicitly known, we impose the solvability condition on the problem (4.26) for </>2 to determine A2 as 

/ w* (L 2 - A ) 4>\ + J w*L 3 w 

Ai — T. . 

J w*w 

Upon using (4.33) for <f>i, L 3 w from (4.23 c), and (L 2 — A ) from (4.23 b), the integrals above are evaluated as 

/ 16tt 2 16A 2 / 4 /- 8 r- ,\ 2tt 4 2 

Ai = - y j c 2 + ^--^r - j C1C0 - — c\ 

V2tt 3y/2ir 4y/2ir n 12^ „ /- /— 

H — c 2 H — c 4 H — c 5 + 8c 6 H — c 7 + 2V2c 8 7r + 2V27rc 9 . 

3 5 3 5 

Finally, upon substituting Co, . . . , eg from (4.24) into this expression, we determine Ai explicitly as 

(4.34, 



2fi 2 t V 72 

The two-term expansion for A given in (4.14) follows from (4.31) and (4.34). The Hopf bifurcation threshold is 
obtained by setting Ai — 0. This occurs precisely as to is increased past ro c , where t c is given by (4.13). ■ 



5 Hot- Spot Patterns in 2-D: Equilibria and Stability 

In this section we construct a if-spot quasi-steady-state solution to (1.1) in an arbitrary 2-D domain with spots 
centered at xi, . . . ,Xk- To leading-order in a = — 1/loge, we then derive a threshold condition on the diffusivity D 
for the stability of the i^-spot quasi-steady-state solution to instabilities that develop on an 0(1) time-scale. 
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As in the analysis of hot-spot patterns in one spatial dimension, we set V = P/A 2 (see (2.2)) into (1.1) to obtain 
A t = s 2 AA- A + VA 3 + a, xefl; d n A = 0, x e dQ , (5.1a) 
r(A 2 V) t = £>V • (A 2 WV) -VA 3 +7 -a, x £ n ; d n V = , xeffi. (5.16) 

We first motivate the e-dependent re-scalings of V and A that are needed for the 2-D case. We suppose that 
D 1, so that V is approximately constant. By integrating the steady-state equation of (5.1 6) over f2 we get 
V = c/ J n A 3 dx, where c is some 0(1) constant. Therefore, if A = (D(s~ p ) in the inner hot-spot region of area 0(e 2 ), 
we obtain J Q A 3 dx = C(e 2 ~ 3p ), so that V — 0(e 3p ~ 2 ). In addition, from the steady-state of (5.1 a), we must have in 
the inner region that A 3 V ~ A, so that -ip + (3p - 2) = -p. This yields p = 2. Therefore, for D > 1, V = 0{e A ) 
globally on ft, while A = 0{e~ 2 ) in the inner region near a hot-spot. Finally, in the outer region we must have 
A ~ a = O(l), so that from (5.1 6), we conclude that DV • (A 2 W) ~ a - 7 = (9(1). Since V = C(e 4 ), this balance 
requires that D = (D(e~ 4 ). Finally, in the core of a hot-spot we conclude that the density P of criminals, given by 
P = VA 2 , is 0(1) ase^O. 

Although this simple scaling analysis correctly identifies the algebraic factors in e, there are more subtle logarithmic 
terms of the form a = —1/ logs that are needed in the construction of the quasi-steady-state hot spot solution. 
The scaling analysis above motivates the introduction of new variables v, u, and V defined by 

V = e A v, A = e- 2 u, D = V/e i . (5.2) 

In terms of (5.2), (5.1) transforms exactly to 

u t = e 2 Au — u + vu 3 + ae 2 , x G £1 ; d n u = , x <G <9f2 , (5.3 a) 

V 

t (u 2 v) t = -^V • (ti 2 Vi>) - e~ 2 vu 3 +7- a, xeQ; d n v = . x e dfl . (5.3 6) 

Owing to the non- uniformity in the behavior as \y\ — > 00 of the solution to this core problem near a spot (see 
below), the construction of a quasi-steady-state if -spot solution for (5.3) is more intricate than that for the Gierer- 
Mcinhardt, Schnakenburg, or Gray-Scott problems analyzed in [39]-[44], [13], [16], and [1]. As such, we will only 
develop a theory that is accurate to leading order in a = — 1/loge, similar to that undertaken in [39]-[44], and 
[13]. This is in contrast to the recent approach in [16] and [1] that used a hybrid asymptotic-numerical method to 
construct quasi-steady-state spot patterns to the Schnakenburg and Gray-Scott systems, respectively, with an error 
that is beyond-all-orders with respect to a. 

In the outer region, away from the spots centered at xi, . . . , xk, we expand u and v as 

u = ae 2 + o(e 2 ) , v ~ h + ah\ + • • • , (5.4) 

where a = — 1/ loge and V = T> /a, where T> = 0(1). From the steady-state of (5.3 b) we obtain that ho is constant, 
and that hi satisfies 

A/ti = - ^"^ , x€ft\{ Xl ,...,x K }; d n h! = 0, xen. (5.5) 
a V 

As shown below, this problem must be augmented by certain singularity conditions that are obtained by matching 
the outer solution for v to certain inner solutions, one in the neighborhood of each spot. 

In the inner region near the j-th spot centered at Xj we introduce the inner variables y, Uj, and Vj by 

y = £- 1 (x-x j ), Vj(y)=v(xj+ey), uj(y) = u( Xj + ey) . (5.6) 
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In terms of these inner variables, and with V = o-~ lr D , the steady-state of (5.3) transforms exactly on y € R 2 to 

AyUj — Uj + VjU 3 - + ae 2 = , (5.7 a) 

4 6 

where a = — 1/loge. We will construct a radially symmetric solution Uj = Uj(p), Vj = vj(p) to this problem, where 

p= \v\- 

The complication in analyzing (5.7) is that uj = 0(1) for |y| = 0(1), whereas Uj = 0(e 2 ) for \y\ ^ 1. Therefore, 
the "diffusivity" u 2 in the operator for v in (5.7 6) ranges from 0(1) when \y\ = 0(1) to 0(e 4 ) when \y\ >• 1. For 
this reason, we cannot simply neglect the second term on the left-hand side of (5.7 6) for all \y\. However, as we show 
below, we can neglect the 0(e 6 ) term on the right-hand side of (5.7 b). 

For \y\ — 0(1), we expand Uj and vj as 

Uj = u j0 + e 2 uji H , vj = Vjo + e 2 Wji +■■■ ■ (5.8) 

Upon substituting this expansion into (5.7), we obtain that Vjo and Vji are constants, and that Uj and Uj\ are 
radially symmetric solutions of 

ApUj - Uj + v j0 u 3 = , ApUji - Uji + 3u 2 j0 v j0 Uji = -a - u 3 Vji , 

on p > with Ujo — > and Uj\ — > a as p — > oo. Here A p g = g" + p _1 <7' for g = .g(p). In terms of the unknown 
constants Vj and Vji, the solutions for u j0 an d Uj\ are 

u jo = v J 1/2w > u ji = a ~TH W ~ 3aw i ' ( 5 - 9 ) 

2v jo 

where w — w(p) and »i = w\(p) are the unique radially symmetric solutions of 

A p w — w + w 3 = , Lw\ = ApWi — wi + 3w 2 wi — w 2 , (5.10) 

with w(p) > 0, w'(0) = 0, and w — > as p — » oo, together with w^(0) = and w\ — > as p — » oo. The expression 
(5.9) for Uji shows that ttji — > a as p — > oo, so that from (5.8) Uj ~ ae 2 when |j/| 3> 1. 

Next, we calculate the far-field behavior, valid for \y\ ^> 1, for the solution Vj to (5.7 6). To do so, we define the 
ball Bs = {y | \y\ < <5}, where 1 <C 5 <C 0(e _1 ). Therefore, this ball is defined in the intermediate matching region 
between the inner and outer scales y and x, respectively Upon integrating (5.7 6) over Bs, and using the divergence 
theorem, we obtain that 

4 r 

2Tru 2 Sv / Ap = s ~ / VjU 3 dy + 0(8 2 ae 6 ) . (5.11) 

U JBs 

Since Uj = 0(1) only for \y\ — 0(1) where Vj ~ w^o + o(l), the integral on the right hand-side of (5.11) can be 
estimated by using Uj <~ v^ 2 w. In contrast, on the left hand-side of (5.11) we use <~ ae 2 on p = S » 1. In this 
way, we obtain that (5.11) becomes 

2na 2 Ss 4 v'j\ p=5 ~ f ^^r") f°° w 3 pdp + 0(6 2 ae 6 ) . (5.12) 
V\A ) 7o' y o/ Jo 

Since <5 <C ©(e^ 1 ), we can neglect the last term on the right hand-side of (5.12), which is equivalent to neglecting 
the right-hand side of (5.7 6) for the inner problem. 
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From (5.12) we obtain that the far- field behavior for |y| » 1 for the solution to (5.7 6) has the form 

Vj ~ v j0 + a (Sj log \y\ + 0(1)) , (5.13 a) 

where Sj is defined by 

b f°° 

Sj = ^ — , b= I w 3 pdp. (5.136) 



Therefore, the appropriate core problem determining the asymptotic shape of the hot-spot profile is to seek a radially 
symmetric solution to 

e 4 a 

A yUj - Uj + Vju) + as 2 = , V„ ■ (tfVyVj) = j^Vju) . (5.14) 

The next step in the construction of the multi hot-spot quasi-steady-state pattern is to match the inner and outer 
solutions for v in order to determine Vjo- We let y = e^ 1 (x — Xj) in (5.13 a) to obtain that the outer solution for v 
must have the singularity behavior 

v ~ Vjo + Sj + a [Sj log | a; — Xj \ + 0(1)] , as x — > Xj , j = 1, . . . , K . (5.15) 

Upon comparing (5.15) with the outer expansion v ~ h + ahi + ■ ■ ■ from (5.4), we conclude that 

h =Vj + Sj, j = l,...,K, (5.16 a) 

and that hi satisfies (5.5) subject to the singularity behaviors 

hi <~ Sjlog \x — Xj\ + 0(1) , as x — > Xj , j = 1, . . . , K . (5.16 6) 

Upon using the divergence theorem, the problem for h\ has a solution only when the solvability condition 



K 



is satisfied, where is the area of fl When this condition is satisfied, the solution for hi can be written as 

K 

hi = -2ir^2 SiG(x; Xi) + hi, (5.17) 

i=l 

where hi is a constant to be determined, and where G(x;Xi) is the Neumann Green's function satisfying 

AG = — !— — 5(x — Xi) , x&fl; d n G = , x e dfl , (5.18 a) 

|S2| 

/ G(x;xi)dx = 0; G ~ — — log \x — Xi\ + 0(1) , as x -> Xj . (5.18 6) 

In summary, the asymptotic matching provides the following algebraic system for determining w,- for j = 1, . . . , if: 

n c _i/2 101(7 — a) 6 

fco=^-o)^o + -=, E^o 1/2 -^^, (5-19) 

V J j = l 

A symmetric if-hot-spot quasi-steady-state solution corresponds to a solution of (5.19) for which Vjo = v for all j. 
From (5.9), this solution is characterized by the fact that the hot-spot profile is, to leading-order, the same for each 
j. The result for such symmetric quasi-equilibria is summarized as follows: 

Principal Result 5.1: For e — > 0, a symmetric K-hot-spot quasi- steady- state solution to (5.3) on the parameter 
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Figure 6. Steady state hot-spot solution of (5.3) in a unit disk. Parameter values are e = 0.05, a = 1,7 = 2,D = 
500e~ 2 ,r = |a;| < 1. (a) The solid line is the steady state solution u(r) of (5.3) computed by solving the associated 
radially symmetric boundary value problem numerically. The dashed line is the asymptotic approximation given by 
(5.20) (b) The solid line is the steady state solution for v(r). Note the "fiat knee" region within the spot center. The 
dashed line is the leading-order asymptotics v ~ vq given by (5.20). 

regime D = e~ a T>q/o- with a = — 1/ loge, is characterized as follows: In the inner region near the j-th hot-spot, where 
y = e (x — Xj), then 

4n 2 b 2 K 2 



u <~ v Q 1 ^ 2 w + 0(e 2 ) , v <~ v a + o(l) . 



vo 



(5.20) 



\fl\ 2 (j~a) 2 ' 

where b = J °° w 3 pdp, and w{p) with p — \y\ is the ground-state solution of (5.10). Alternatively, in the outer region 
where \x — Xj\ >• O(e) for j = 1, . . . , K, then 

u <~ ae 2 , v <~ h + ahi + • • • . (5-21) 

Here h\ is given in (5.17) in terms of the Neumann Green's function, and h is a constant given by 

b An 2 b 2 K 2 |fi|(7-a) 



h = v 



(5.22) 



a 2 V ^/vE \fl\ 2 (j-a) 2 2ita 2 KV ' 
For a symmetric hot-spot pattern, the source strength Sj in (5.17) is the same for each j, and is given in terms of 
vq by Sj = So = bj (a 2 2?ovW) where b = Jg 00 w i pdp. 

To illustrate Principal Result 5.1, we let be the unit disk centered at the origin, and we choose a = 1, 7 = 2, e = 
0.05, and D — 500e -2 . We consider a single hot-spot at the center of the disk. Upon numerically computing the 
ground-state solution w satisfying (5.10) we obtain that 



w(Q) w 2.2062, 
The asymptotic result (5.20) then yields 



w 3 dy w 15.1097. 



v 



(15.1097) 2 

(7 — a) it 2 



23.13184; u(0) ~ w(0)vo 



1/2 



0.4587. 
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Alternatively, from the full numerical solution of the radially symmetric steady-state solution of (5.3), we compute 
that v(0) w 23.017 and u(0) w 0.455. The error is about 0.5% for v(0) and about 0.8% for u(0). A comparison of 
asymptotic and numerical results is shown in Fig. 6. 

As a remark, the general shape of the function T{v) defined in (5.19) also shows that there can be asymmetric 
-ftf-spot quasi-equilibria corresponding to spots of two distinct heights. Similar asymmetric patterns in 2-D have been 
constructed for the Gierer-Meihnardt and Gray-Scott systems in [43]) and ([44]. Let K s and Kb be no n- negative 
integers denoting the number of small and large spots, respectively, with K = K s + Kb- Then, from (5.19), a K-spot 
asymmetric pattern is constructed by determining two distinct values Vo s and Vob satisfying 

F{Vm) = F{v 0r ) , V ah < Wmin < v 0s , -= + -= = — , (5.23) 

y'VQs V«06 27T0 

where w m ; n = 3 (c/2) 3 ^ 2 and c is defined in (5.19). The spatial profile of the hot-spot of large and small amplitude is 
given by u <~ Vq^ 2 w and u <~ Vq^ 2 w, respectively. 

We will not investigate the solvability with respect to T> of the algebraic system (5.23) governing asymmetric 
spot patterns. Instead, in the next subsection we will study the stability properties of the symmetric AT-hot-spot 
quasi-steady-state solution given in Principal Result 5.1. 



5.1 The Stability of Hot-Spot Quasi-Steady-State Patterns 

We linearize (5.3) around the quasi-steady-state if-hot-spot pattern to obtain the eigenvalue problem 



s 2 A(p -<t>+ 3u 2 v 4> + u 3 tp = \<f> , 
1 



x eft: 



d n <f> = . 



(5.24 a) 



DV ■ (u 2 VTp + 2ucl)Vv) ~ — (3u 2 v(t> + u 3 ip) =T\(u 2 ip + 2uv(t>) , xefl; d n ip = , x e dfl . (5.24 6) 

In our stability analysis we will consider the range of D where D = e~ 4 T> /<J and Do = 0(1). It is on this parameter 
range of D that a stability threshold occurs. Our main stability result is as follows: 

Principal Result 5.2: Consider the K-spot quasi- steady- state solution of (5.3) as constructed in Principal Result 
5.1 for e >C 1. Assume that r = tq/e 2 where tq — 0(1). Then, for e — > 0, and to leading order in a = — 1/loge, the 
stability on an Oil) time-scale of a K-hot-spot quasi- steady- state solution with K > 2 is determined by the spectrum 
of the two distinct NLEP's 



A y $ - $ + 3w 2 $ - KiW 3 



J R2 w 2 <S>dy + 2 ^/ Ra dy 
Jr2 w 3 dy 3 0/> ' 



Jr 2 w 2 dy 



A*, 



ye 



with $ — > as \y\ — > oo. Here m for i = 1, 2 are defined by 



Ki = 3(l + r A^) 1 , K 2 = 3(l + r A^ + 



2V 0c 



in terms of the constants f3 and T>q c defined by 

K J R2 w 2 dy 



|fi|(7" 



V 0c 



\Cl\ 3 ( 7 -aY 



v2 ■ 



(5.25 a) 



(5.25 6) 



(5.25 c) 



AiroPK 3 (/ R2 w 3 rfy)' 

For a one-hot-spot solution, for which K = 1, t/iere is oriiy a single NLEP with m in (5.25 a) replaced by k\. For 
t <C 1, w/wc/i is equivalent to t <C ©(e -2 ), we conclude that a K-spot pattern with K > 2 is stable when T> < T> 0c 
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and is unstable when T> n > 2? 0c - In terms of the unsealed D, this yields the stability threshold 

D=(—}v 0c , for K>2, a = -1/ logs. (5.26) 



a 

For K = 1 and To <C 1, a single hot-spot is stable for all Vq independent of e. 

To show that T>q c is the stability threshold of T>o for K > 2, we let tq — > in (5.25) to obtain the limiting NLEP's 

A„$ - $ + 3w 2 $ - ct 3 W f = A$ , i/el 2 : $^0 as Id oo , (5.27) 

J R2 w 3 dj/ 

where k can assume either of the two values 

Kio = 3, K20 = 3^1 + ^-^ . (5.28) 

For (5.27), the rigorous result in Theorem 3 of [45] establishes the existence of an eigenvalue with Re(A) > whenever 
k < 2. Thus, since K20 < 2 when T>o > T>oc, we conclude that a if-hot-spot quasi-equilibria with K > 2 is unstable 
when 2?o > D Oc- In addition, the result in Theorem 1 of [45] (see the remark following Theorem 1) proves that all 
eigenvalues satisfy Re(A) < when 2 < k < 6. Since kiq — 3 and 2 < K20 < 3 for all Vq < T>q Ci it follows that the 
NLEP for a if-hot-spot quasi-steady-state solution with K > 2 does not have unstable eigenvalues when 2?o < T>q c . 
For K = 1, which corresponds to a one- hot-spot solution, the only choice for /c is k — Kio = 3, and so we have 
stability for any T> > independent of e. 

We make two remarks. Firstly, we anticipate that a single hot-spot solution where K = 1 will be stable provided 
that Do is not exponentially large in 1/e. The analysis to determine this stability threshold in the near-shadow limit 
should be similar to that done for the Gierer-Meinhardt model in [14]. Secondly, it is an open question to analyze 
(5.25) for tq — 0(1) to determine if there are any Hopf bifurcations. The analysis of this problem in the 2-D context 
is much more difficult than in 1-D since the identity Lqw 2 = 3u> 2 for the local operator Lq^ = Aj,<£> — $ + 3w 2 $ in M. 1 
no longer holds in R 2 . Moreover, since this identity does not hold in the 2-D case, we cannot determine A explicitly 
when r < 1 as was done for the 1-D case in Lemma 3.2. 

The stability threshold in Principal Result 5.2 follows once we derive the NLEP (5.25). The derivation of this 
NLEP is done in several distinct steps. 

We first consider the inner region near the j-th spot and we introduce the new variables y, <&j(y), and ^ j{y) by 

V = s~ 1 (x-x j ), $>j(y) = 4>{xj + sy) , V j (y) = i>(x j +ey). 
Then, with D = e^ 4 V /a and a = -1/loge, (5.24) on y e R 2 becomes 

- '!'.. + *j + " : ; *j - -\<!>, , (5-29 a) 

V y ■ (tfVyVj + 2u j $ j V y v j ) - (3u 2 jVj + u] 9 j) = T -^- (up j + 2u j v j $ j ) . (5.29 6) 

Since uj = O(l) and V y Vj = o(l) when \y\ = 0(1), we obtain from (5.29 b) that ^ = V j0 + O(e 2 ) when \y\ = 0(1), 
where ^S/jo is an unknown constant. Then, with Uj <~ w/y/vo an d Vj <~ v for \y\ = 0(1), we obtain from (5.29 a) that 
<&j ~ <F,o + o(l), where $^0 satisfies 

Aj / $ j0 -^o + 3w 2 $ J o + V 3/2 w 3 * j0 = A*jo, i/et 2 ; $ i0 4 as |y|->-oo, (5.30) 

where w is the radially symmetric ground-state solution satisfying (5.10). 
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The next step in the analysis is to determine the constant >3/jo m (5.30). This is done by first determining the 
far- field behavior as \y\ — > oo of the solution to (5.29 b). We define the ball B 5 = {y \ \y\ < 5}, where 1 < 5 < C(£ _1 ), 
and we integrate (5.29 b) over Bg to obtain 

(2wup' j \ p=s + Anu^jv'^s) S = £ -f f (* ,u) + 3u 2 ^) rfy + 1^ [ + 2t W $ j ) . (5.31) 

We now estimate the terms in (5.31) for e <C 1. 

Since the dominant contribution to the integrals on the right hand-side of (5.31) arises from the region where 
\y\ = 0(1), we can asymptotically estimate these integrals by using Uj ~ w/^/vo, Vj <~ t>o, and <~ "J^o- For the 
left hand-side of (5.31) we use Uj ~ ae 2 on p = |j/| = <5 >• 1 to get 



e 4 a 



(2TTe 4 a 2 ^' J \ p=s +47re 2 a^$ 3 | p=(5 ) ^ ^ ^ 



3/2 



/ w 3 dy + 3 w 2 $ j0 dy 

JR 2 JR 



+ ■ 



rAe 6 



V JR2 JR2 



(5.32) 



Next, we estimate the second term on the left hand-side of (5.32). For \y\ » 1, we obtain from the outer limit of 
(5.29 a) that + e 3 a 6 ^j ~ so that with ^ - + o(l), we get 



4\ 



e 6 a 3 



on p = |y| = (5 > 1 



(5.33) 



3 1 + A 

In addition, from (5.13 a), we estimate that v'j\ p =s ~ crSo/5, where So is defined in Principal Result 5.1. Substituting 
this estimate together with (5.33) into (5.32) we obtain 



. . . 2a 2 e 4 aSo 



2irT> a 2 



"*£ f 

yl' 2 



w A dy + 3 / w 2 $ j0 dy 



r- 



/ w 2 dy + 2 v ^ [ w$> lQ dy 

V JR2 JR2 



(5.34) 



2-KV a 2 

From (5.34), and under the assumption that r = £~ 2 t , where r = 0(1), we conclude that tyj has the far-field 
behavior 



j - j0 + <x log + 0(1)) , for |y| » 1 



where Bj for j = 1, . . . , K is defined by 



Si 



1 



27rX> a 2 



jo 



3/2 



w 3 dy + 3 / w 2 ®^ dy 



t A 



/ w 2 dy + 2^ a f w<S> lQ dy 

V JR2 JR2 



(5.35) 



(5.36) 



r2 J R 2 2irT> a 2 

Upon writing (5.35) in terms of the outer variable (x — Xj) = ey, we obtain that the matching condition for the outer 
solution ip is 

tp ~ Bj + $ j0 + a(Bjlog\x - Xj\ + 0(1)) , as x -> Xj , j = 1, . . . (5.37) 

Next, we consider the outer region for ip where \x — Xj\ = 0(1) for j = 1,...,K. In (5.24 b) we use u <~ ae 2 , 
Vjj = 0(cr), and = 0(e 6 ) from (5.33) to estimate that 

• (oVV^ + 0(ae 8 )) - e~ 2 (a 3 s 6 ^ + 0(ae 10 )) = rA (a 2 eV + 0(ae» . 
Hence, when r = £~ 2 t we obtain that T> Aip ~ To\s 2 aip. In order to match with (5.37) we must expand the outer 
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solution for tp as 

V> = Vo + <n/>i + ■ ■ ■ • (5.38) 

We then obtain that ip is a constant, given by 

i, = ^ j0 + B j , j = l,...,K, (5.39) 

and that satisfies 

AV>i=0, x Gft\{xi,...,x K }; d„ip! = , xeffl, (5.40 a) 

V>i ~ Bj log\x - Xj\ + (9(1) , as x ^ Xj , j = 1,...,K . (5.40 b) 

The solvability condition for (5.40) is that J2f=i Bj = 0. Upon summing (5.39) from j = 1, . . . , K, we obtain that 

1 K 

B > ■ = -*><>+ k^Vjo, (5-41) 

where -B^ is given in (5.36). 

The final step in the derivation of the NLEP is to solve (5.36) and (5.41) for and substitute the resulting 

expression into (5.30). To do so, We introduce the vectors B = [B\, . . . , Bk) T and = (^oi, • • • , ^ok) T , and we 
write the system (5.36) and (5.41) in matrix form as 

B = ci^o - C2T0 + c 3 *o - C4.F1 , B = -(J-£)tf , (5.42) 

where the constants Ci for i = 1, . . . , 4 are defined by 

ci = — ^ f w 3 dy, c 2 = 3ci«o /2 , c 3 = - J"° A f w 2 dy , c 4 = 2uo /2 c 3 , (5.43 a) 

27rV a a 2 v 3 /2 V 2irV a 2 v J R2 

and the vectors Fq and J"i are defined by 

r - /r 2 w2 ®° d y Tl - A 2 d y (5 4 36 ) 

° J R 2W 3 dy ' 1 J R 2W 2 dy 

Here l»o = ($01, • • • j ®ok) T ■ In (5.42), I is the K x K identity matrix and the matrix £ is defined by £ — K~ 1 ee T , 
where e = (1, . . . , 1) T . 

By solving (5.42) for ^0, and substituting the resulting expression into (5.30), we obtain the vector NLEP 



A„$ - $ + 3u/$o + w s M [ T + -^Ti J = A$ , (5-44) 



where the matrix M. is defined by 



(l + c 1 +c 3 ) T 1 



-1 



I-—E . (5.45) 
c 2 c 2 

The matrix M^ 1 is a rank-one update of a scalar multiple of the identity matrix. As such, its spectrum Muj = kuj 
can readily be calculated as 

Kl = ~ — 3/2' = (1,...,1) T ; k 2 = ... = k k = 372, wje = 0, j = 2,...,K. 

(ci + c 3 )v Q ' (1 + ci + c 3 )w 

Notice that the eigenvectors corresponding to the matrix eigenvalues kj for j = 2, . . . , K span the if — 1 dimensional 
subspace perpendicular to e = (1, . . . , 1) T . As such, the competition instability modes correspond to j — 2, . . . , K, 
whereas the synchronous instability mode corresponds to k,\ with eigenvector lo\ = (1, . . . , 1) T . 
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Figure 7. Numerically computed bifurcation diagram of A(0) vs. 7. The parameter values are a = l,e = 0.05, a; G 
[0, 1], and D = 2. A localized hot-spot appears for large values of A(0). The asymptotics A(0) <~ 2< -g~"' > (sec (2.19)) 
are shown by a dotted line. The constant steady state A ~ 7 is indicated by a solid straight line line. Turing patterns 
are born from the spatially uniform steady state as a result of a Turing bifurcation at 7 ~ 3a/2 = 1.5. The weakly 
nonlinear regime is indicated by a dashed parabola coming out of the bifurcation point. Inserts shows the change in 
the shape of the profile A(x) along the bifurcation curve. 

Then, we use the explicit formulae for cj in (5.43 a) and for vq in (5.20) to write Ki for i = 1, 2 as in (5.25 b). In 
addition, the ratio C4/C2 in (5.44) can be calculated using (5.43 a) and (5.20) to get C4/C2 = 2toA/3/3, where ft is 
defined in (5.25 c). Finally, by diagonalizing the vector NLEP (5.44) by using the matrix decomposition of M., and 
by recalling the definition of Ti for i = 0, 1 in (5.43 b), we obtain the NLEP (5.25) of Principal Result 5.2. This 
completes the derivation of Principal Result 5.2 ■. 

5 Discussion 

We have studied localized hot-spot solutions of (1.1) in one and two spatial dimensions in the regime e 2 -C 1 with 
D ^> 1. In this large D limit, steady-state multi hot-spot solutions have been constructed and their stability properties 
investigated with respect to D and r from the analysis of certain nonlocal eigenvalue problems. An open problem is to 
characterize the dynamics of multi hot-spot patterns by reducing (1.1) to a finite dimensional dynamical system for 
the locations of the hot-spots in a quasi-steady-state pattern as was done for various two-component reaction-diffusion 
models without drift terms in [6, 7, 9, 33, 16, 1]. 

We now remark on how the localized states constructed in this paper are related to the weakly nonlinear Turing 
patterns studied in [29, 30, 31]. In contrast to the theory developed in [30, 31], our parameter values for (1.1) arc 
not restricted to lie close to the Turing bifurcation point. In the limit e — > 0, the spatially homogeneous steady-state 
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solution for (1.1) is A e = 7 and P e = (7 — a)/-f. For e — > 0, it is linearly unstable when (see equation (2.8) of [30]) 

3 

7 > —a , as e — > , 

and a spatially heterogeneous solution bifurcates off the homogeneous steady state at 7 ~ |a as e — > 0. In contrast, a 
localized hot-spot exists for the wider parameter range 7 > a. In Fig. 7 we plot the numerically computed bifurcation 
diagram for ^4(0) vs. 7 on a one-dimensional interval, with other parameters as indicated in the figure caption. The 
Turing bifurcation at 7 = |a is subcritical when e <C 1 (cf. [30]), which is consistent with the stability of the constant 
state when 7 < |a. However, the bifurcation curve quickly turns around as it enters the localized regime. This is 
consistent with the existence of localized states when 7 > a. 

The dispersion relation obtained by linearizing (1.1) about the spatially uniform state A e , P e is calculated as 

tA 2 + A (A e + Dm 2 + t(1 - P e ) + re 2 m 2 ) + A e (l + em 2 ) - 3DP e m 2 + Dm 2 {\ + s 2 m 2 ) = , (5.1) 

where A e — 7 and P e = (7 — a) /j. For 7 > |a, it is readily shown from this relation that the edges of the Turing 
instability band miower < m < m uppcr satisfy 

mi owor ~ L>- 1 / 2 7(2 7 - 3a)- 1 / 2 ; m uppcr ~ e^" 1 / 2 ^ - 3a) 1 / 2 , as e^O. 

The most unstable mode (i.e. the one which grows the fastest, and therefore the one most commonly observed) is 
obtained by setting dX/dm = in (5.1). For e — > 0, this gives the maximum growth rate 

Adominant ^ 3-P e 1 2s TTl , 

together with the most unstable mode 

m dominant ~ e-^D- 1 '^- 1 ' 2 [(7 - a)(3 7 2 + 2r(2 7 - 3a)] 1/4 , 

which is consistent with [29] in the limit e — >• (see also equation (2.9) of [30]). 

Correspondingly, this implies that for an initial condition consisting of a random perturbation of the spatially 
uniform steady- state, the preferred pattern has a characteristic half-length during ~ 7r/mdominant, where 

Zturing - e^ 2 D^W 12 [(7 - «)(37 2 + 2r(2 7 - 3a))] " 1/4 n . 

In contrast, for localized structures the characteristic length I between hot-spots to ensure stability of a multi hot-spot 
pattern, as obtained from (1.4), is that I > l c where 

V^D 1/4 £ 1/2 a 1/2 (7-a)- 3 / 4 . (5.2) 

Although during and l c are not related, they are both of the same asymptotic order 0(D 1 ^ 4 e 1 ^ 2 ), which implies 
that the number of stable localized hot-spots corresponds roughly to the most unstable Turing mode. In fact, for a 
large parameter range, the inequality l c < during holds. For example, when r = l,a = 1,7 = 2 we calculate that 
lc/ during = 0.7 < 1. This was already observed empirically in Fig. 6 of [29]. 

Formula (1.4) provides an upper bound K c on the number K of stable hot-spots in the regime where D 1. 
On the other hand, numerical evidence shows the existence of a lower bound that occurs when D is sufficiently 
small. In this regime an instability occurs when there are too few hot-spots. In fact, if the hot-spot inter-distance 
exceeds some critical length, then a hot-spot insertion phenomena is observed (see Fig. 8). From Fig. 8, it appears 
that hot-spot insertion takes place every time that D is quartered. A similar insertion phenomenon occurs in other 
reaction-diffusion systems such as a chemotaxis model [26], the Gray-Scott model [27], the ferrocyanide-iodide-sulfitc 
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Figure 8. Hot-spot insertion phenomenon. Numerical solution of (1.1) with e — 0.02, a=l,7 = 2,D = l/(l + O.Oli) 
with x G (0, 2). Initial conditions consist of two boundary hot-spots. Snapshots of P(x) are shown for values of D as 
indicated. Hot-spot insertion takes place every time that D is quartered. 



system [20], the Brusselator [17], the Schnakenburg model [16], and a model of droplet breakup [21]. The detailed 
analysis of hot-spot insertion phenomena in (1.1) will be considered in future work. 
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